/************************************************************************
 * $Id: astro.c,v 1.8 2004/12/10 14:08:04 thamer Exp $
 *
 * ------------
 * Description:
 * ------------
 *  Copyright (c) 2003, Arabeyes, Thamer Mahmoud
 *
 *  A full featured Muslim Prayer Times calculator
 *
 *  Most of the astronomical values and formulas used in this file are
 *  based on the Jean Meeus subset of the VSOP87 planetary theory.
 *
 * -----------------
 * Revision Details:    (Updated by Revision Control System)
 * -----------------
 *  $Date: 2004/12/10 14:08:04 $
 *  $Author: thamer $
 *  $Revision: 1.8 $
 *  $Source: /home/arabeyes/cvs/projects/itl/libs/prayertime/src/astro.c,v $
 *
 * (www.arabeyes.org - under LGPL license - see COPYING file)
 ************************************************************************/

#include "astro.h"


typedef struct
{
    double dec;
    double ra;
    double sidtime;
    double dra;
    double rsum;
        
} AstroDay ;


static void computeAstroDay(double JD, AstroDay* astroday);
static void computeTopAstro(const Location* loc, const Astro* astro, Astro* topAstro);

double getRefraction(const Location* loc, double sunAlt)
{
    double part1, part2;

    part1 = (loc->pressure/1010.0) * (283/(273 + loc->temperature));
    part2 = 1 / tan(DEG_TO_RAD(sunAlt + (7.31/(sunAlt + 4.4)))) + 0.0013515;

    return (part1 * part2) / 60.0;
}

double getJulianDay(const Date* date, double gmt)
{
    double jdB=0, jdY, jdM, JD;
    
    jdY=date->year;
    jdM=date->month;

    if (date->month <= 2) {
        jdY--;
        jdM+=12;
    }

    if (date->year < 1) 
        jdY++;

    if ((date->year > 1582) || ((date->year == 1582) &&
                                ((date->month > 10) ||
                                 ((date->month == 10) && (date->day >= 4)))))
        jdB = 2 - floor(jdY/100.0) + floor((jdY/100.0)/4.0);

    JD = floor(365.25 * (jdY + 4716.0)) + floor(30.6001 * ( jdM + 1))
        + (date->day + (-gmt)/24.0)  + jdB - 1524.5 ;

    return JD;

}

const double L0[64][3]={
    {175347046, 0, 0},
    {3341656, 4.6692568, 6283.07585},
    {34894, 4.6261, 12566.1517},
    {3497, 2.7441, 5753.3849},
    {3418, 2.8289, 3.5231},
    {3136, 3.6277, 77713.7715},
    {2676, 4.4181, 7860.4194},
    {2343, 6.1352, 3930.2097},
    {1324, 0.7425, 11506.7698},
    {1273, 2.0371, 529.691},
    {1199, 1.1096, 1577.3435},
    {990, 5.233, 5884.927},
    {902, 2.045, 26.298},
    {857, 3.508, 398.149},
    {780, 1.179, 5223.694},
    {753, 2.533, 5507.553},
    {505, 4.583, 18849.228},
    {492, 4.205, 775.523},
    {357, 2.92, 0.067},
    {317, 5.849, 11790.629},
    {284, 1.899, 796.298},
    {271, 0.315, 10977.079},
    {243, 0.345, 5486.778},
    {206, 4.806, 2544.314},
    {205, 1.869, 5573.143},
    {202, 2.4458, 6069.777},
    {156, 0.833, 213.299},
    {132, 3.411, 2942.463},
    {126, 1.083, 20.775},
    {115, 0.645, 0.98},
    {103, 0.636, 4694.003},
    {102, 0.976, 15720.839},
    {102, 4.267, 7.114},
    {99, 6.21, 2146.17},
    {98, 0.68, 155.42},
    {86, 5.98, 161000.69},
    {85, 1.3, 6275.96},
    {85, 3.67, 71430.7},
    {80, 1.81, 17260.15},
    {79, 3.04, 12036.46},
    {71, 1.76, 5088.63},
    {74, 3.5, 3154.69},
    {74, 4.68, 801.82},
    {70, 0.83, 9437.76},
    {62, 3.98, 8827.39},
    {61, 1.82, 7084.9},
    {57, 2.78, 6286.6},
    {56, 4.39, 14143.5},
    {56, 3.47, 6279.55},
    {52, 0.19, 12139.55},
    {52, 1.33, 1748.02},
    {51, 0.28, 5856.48},
    {49, 0.49, 1194.45},
    {41, 5.37, 8429.24},
    {41, 2.4, 19651.05},
    {39, 6.17, 10447.39},
    {37, 6.04, 10213.29},
    {37, 2.57, 1059.38},
    {36, 1.71, 2352.87},
    {36, 1.78, 6812.77},
    {33, 0.59, 17789.85},
    {30, 0.44, 83996.85},
    {30, 2.74, 1349.87},
    {25, 3.16, 4690.48}
};

const double L1[][3]={
    {628331966747.0, 0, 0},
    {206059, 2.678235, 6283.07585},
    {4303, 2.6351, 12566.1517},
    {425, 1.59, 3.523},
    {119, 5.796, 26.298},
    {109, 2.966, 1577.344},
    {93, 2.59, 18849.23},
    {72, 1.14, 529.69},
    {68, 1.87, 398.15},
    {67, 4.41, 5507.55},
    {59, 2.89, 5223.69},
    {56, 2.17, 155.42},
    {45, 0.4, 796.3},
    {36, 0.47, 775.52},
    {29, 2.65, 7.11},
    {21, 5.34, 0.98},
    {19, 1.85, 5486.78},
    {19, 4.97, 213.3},
    {17, 2.99, 6275.96},
    {16, 0.03, 2544.31},
    {16, 1.43, 2146.17},
    {15, 1.21, 10977.08},
    {12, 2.83, 1748.02},
    {12, 3.26, 5088.63},
    {12, 5.27, 1194.45},
    {12, 2.08, 4694},
    {11, 0.77, 553.57},
    {10, 1.3, 3286.6},
    {10, 4.24, 1349.87},
    {9, 2.7, 242.73},
    {9, 5.64, 951.72},
    {8, 5.3, 2352.87},
    {6, 2.65, 9437.76},
    {6, 4.67, 4690.48}
};

const double L2[][3]={
    {52919, 0, 0},
    {8720, 1.0721, 6283.0758},
    {309, 0.867, 12566.152},
    {27, 0.05, 3.52},
    {16, 5.19, 26.3},
    {16, 3.68, 155.42},
    {10, 0.76, 18849.23},
    {9, 2.06, 77713.77},
    {7, 0.83, 775.52},
    {5, 4.66, 1577.34},
    {4, 1.03, 7.11},
    {4, 3.44, 5573.14},
    {3, 5.14, 796.3},
    {3, 6.05, 5507.55},
    {3, 1.19, 242.73},
    {3, 6.12, 529.69},
    {3, 0.31, 398.15},
    {3, 2.28, 553.57},
    {2, 4.38, 5223.69},
    {2, 3.75, 0.98}
};

const double L3[][3]={
    {289, 5.844, 6283.076},
    {35, 0, 0},
    {17, 5.49, 12566.15},
    {3, 5.2, 155.42},
    {1, 4.72, 3.52},
    {1, 5.3, 18849.23},
    {1, 5.97, 242.73}   
};

const double L4[][3]={
    {114.0, 3.142, 0.0},
    {8.0, 4.13, 6283.08},
    {1.0, 3.84, 12566.15}
};

const double L5[][3]={
    {1, 3.14, 0}
};

const double B0[][3]={

    {280, 3.199, 84334.662},
    {102, 5.422, 5507.553},
    {80, 3.88, 5223.69},
    {44, 3.7, 2352.87},
    {32, 4, 1577.34}
};

const double B1[][3]={

    {9, 3.9, 5507.55},
    {6, 1.73, 5223.69}
};

const double R0[][3]={
    {100013989, 0, 0},
    {1670700, 3.0984635, 6283.07585},
    {13956, 3.05525, 12566.1517},
    {3084, 5.1985, 77713.7715},
    {1628, 1.1739, 5753.3849},
    {1576, 2.8469, 7860.4194},
    {925, 5.453, 11506.77},
    {542, 4.564, 3930.21},
    {472, 3.661, 5884.927},
    {346, 0.964, 5507.553},
    {329, 5.9, 5223.694},
    {307, 0.299, 5573.143},
    {243, 4.273, 11790.629},
    {212, 5.847, 1577.344},
    {186, 5.022, 10977.079},
    {175, 3.012, 18849.228},
    {110, 5.055, 5486.778},
    {98, 0.89, 6069.78},
    {86, 5.69, 15720.84},
    {86, 1.27, 161000.69},
    {85, 0.27, 17260.15},
    {63, 0.92, 529.69},
    {57, 2.01, 83996.85},
    {56, 5.24, 71430.7},
    {49, 3.25, 2544.31},
    {47, 2.58, 775.52},
    {45, 5.54, 9437.76},
    {43, 6.01, 6275.96},
    {39, 5.36, 4694},
    {38, 2.39, 8827.39},
    {37, 0.83, 19651.05},
    {37, 4.9, 12139.55},
    {36, 1.67, 12036.46},
    {35, 1.84, 2942.46},
    {33, 0.24, 7084.9},
    {32, 0.18, 5088.63},
    {32, 1.78, 398.15},
    {28, 1.21, 6286.6},
    {28, 1.9, 6279.55},
    {26, 4.59, 10447.39}
};

const double R1[][3]={

    {103019, 1.10749, 6283.07585},
    {1721, 1.0644, 12566.1517},
    {702, 3.142, 0},
    {32, 1.02, 18849.23},
    {31, 2.84, 5507.55},
    {25, 1.32, 5223.69},
    {18, 1.42, 1577.34},
    {10, 5.91, 10977.08},
    {9, 1.42, 6275.96},
    {9, 0.27, 5486.78}
};

const double R2[][3]={

    {4359, 5.7846, 6283.0758},
    {124, 5.579, 12566.152},
    {12, 3.14, 0},
    {9, 3.63, 77713.77},
    {6, 1.87, 5573.14},
    {3, 5.47, 18849}

};

const double R3[][3]={
    {145, 4.273, 6283.076},
    {7, 3.92, 12566.15}
};

const double R4[][3]={
    {4, 2.56, 6283.08}
};

const double PE[][4]={
    {-171996, -174.2, 92025, 8.9},
    {-13187, -1.6, 5736, -3.1},
    {-2274, -0.2, 977, -0.5},
    {2062, 0.2, -895, 0.5},
    {1426, -3.4, 54, -0.1},
    {712, 0.1, -7, 0},
    {-517, 1.2, 224, -0.6},
    {-386, -0.4, 200, 0},
    {-301, 0, 129, -0.1},
    {217, -0.5, -95, 0.3},
    {-158, 0, 0, 0},
    {129, 0.1, -70, 0},
    {123, 0, -53, 0},
    {63, 0, 0, 0},
    {63, 0.1, -33, 0},
    {-59, 0, 26, 0},
    {-58, -0.1, 32, 0},
    {-51, 0, 27, 0},
    {48, 0, 0, 0},
    {46, 0, -24, 0},
    {-38, 0, 16, 0},
    {-31, 0, 13, 0},
    {29, 0, 0, 0},
    {29, 0, -12, 0},
    {26, 0, 0, 0},
    {-22, 0, 0, 0},
    {21, 0, -10, 0},
    {17, -0.1, 0, 0},
    {16, 0, -8, 0},
    {-16, 0.1, 7, 0},
    {-15, 0, 9, 0},
    {-13, 0, 7, 0},
    {-12, 0, 6, 0},
    {11, 0, 0, 0},
    {-10, 0, 5, 0},
    {-8, 0, 3, 0},
    {7, 0, -3, 0},
    {-7, 0, 0, 0},
    {-7, 0, 3, 0},
    {-7, 0, 3, 0},
    {6, 0, 0, 0},
    {6, 0, -3, 0},
    {6, 0, -3, 0},
    {-6, 0, 3, 0},
    {-6, 0, 3, 0},
    {5, 0, 0, 0},
    {-5, 0, 3, 0},
    {-5, 0, 3, 0},
    {-5, 0, 3, 0},
    {4, 0, 0, 0},
    {4, 0, 0, 0},
    {4, 0, 0, 0},
    {-4, 0, 0, 0},
    {-4, 0, 0, 0},
    {-4, 0, 0, 0},
    {3, 0, 0, 0},
    {-3, 0, 0, 0},
    {-3, 0, 0, 0},
    {-3, 0, 0, 0},
    {-3, 0, 0, 0},
    {-3, 0, 0, 0},
    {-3, 0, 0, 0},
    {-3, 0, 0, 0}
};

const int SINCOEFF[][5]={
    {0, 0, 0, 0, 1},
    {-2, 0, 0, 2, 2},
    {0, 0, 0, 2, 2},
    {0, 0, 0, 0, 2},
    {0, 1, 0, 0, 0},
    {0, 0, 1, 0, 0},
    {-2, 1, 0, 2, 2},
    {0, 0, 0, 2, 1},
    {0, 0, 1, 2, 2},
    {-2, -1, 0, 2, 2},
    {-2, 0, 1, 0, 0},
    {-2, 0, 0, 2, 1},
    {0, 0, -1, 2, 2},
    {2, 0, 0, 0, 0},
    {0, 0, 1, 0, 1},
    {2, 0, -1, 2, 2},
    {0, 0, -1, 0, 1},
    {0, 0, 1, 2, 1},
    {-2, 0, 2, 0, 0},
    {0, 0, -2, 2, 1},
    {2, 0, 0, 2, 2},
    {0, 0, 2, 2, 2},
    {0, 0, 2, 0, 0},
    {-2, 0, 1, 2, 2},
    {0, 0, 0, 2, 0},
    {-2, 0, 0, 2, 0},
    {0, 0, -1, 2, 1},
    {0, 2, 0, 0, 0},
    {2, 0, -1, 0, 1},
    {-2, 2, 0, 2, 2},
    {0, 1, 0, 0, 1},
    {-2, 0, 1, 0, 1},
    {0, -1, 0, 0, 1},
    {0, 0, 2, -2, 0},
    {2, 0, -1, 2, 1},
    {2, 0, 1, 2, 2},
    {0, 1, 0, 2, 2},
    {-2, 1, 1, 0, 0},
    {0, -1, 0, 2, 2},
    {2, 0, 0, 2, 1},
    {2, 0, 1, 0, 0},
    {-2, 0, 2, 2, 2},
    {-2, 0, 1, 2, 1},
    {2, 0, -2, 0, 1},
    {2, 0, 0, 0, 1},
    {0, -1, 1, 0, 0},
    {-2, -1, 0, 2, 1},
    {-2, 0, 0, 0, 1},
    {0, 0, 2, 2, 1},
    {-2, 0, 2, 0, 1},
    {-2, 1, 0, 2, 1},
    {0, 0, 1, -2, 0},
    {-1, 0, 1, 0, 0},
    {-2, 1, 0, 0, 0},
    {1, 0, 0, 0, 0},
    {0, 0, 1, 2, 0},
    {0, 0, -2, 2, 2},
    {-1, -1, 1, 0, 0},
    {0, 1, 1, 0, 0},
    {0, -1, 1, 2, 2},
    {2, -1, -1, 2, 2},
    {0, 0, 3, 2, 2},
    {2, -1, 0, 2, 2}
};

void getAstroValuesByDay(double julianDay, const Location* loc, Astro* astro, 
                         Astro* topAstro)
{
    AstroDay ad; 

    if (astro->jd == julianDay-1)
    {
        astro->ra[0] = astro->ra[1];
        astro->ra[1] = astro->ra[2];
        astro->dec[0] = astro->dec[1];
        astro->dec[1] = astro->dec[2];
        astro->sid[0] = astro->sid[1];
        astro->sid[1] = astro->sid[2];
        astro->dra[0] = astro->dra[1];
        astro->dra[1] = astro->dra[2];
        astro->rsum[0] = astro->rsum[1];
        astro->rsum[1] = astro->rsum[2];
        computeAstroDay(julianDay+1, &ad);
        astro->ra[2] = ad.ra;
        astro->dec[2] = ad.dec;
        astro->sid[2] = ad.sidtime;
        astro->dra[2] = ad.dra;
        astro->rsum[2] = ad.rsum;
    }
    else if (astro->jd == julianDay + 1)
    {

        astro->ra[2] = astro->ra[1];
        astro->ra[1] = astro->ra[0];
        astro->dec[2] = astro->dec[1];
        astro->dec[1] = astro->dec[0];
        astro->sid[2] = astro->sid[1];
        astro->sid[1] = astro->sid[0];
        astro->dra[2] = astro->dra[1];
        astro->dra[1] = astro->dra[0];
        astro->rsum[2] = astro->rsum[1];
        astro->rsum[1] = astro->rsum[0];
        computeAstroDay(julianDay-1, &ad);
        astro->ra[0] = ad.ra;
        astro->dec[0] = ad.dec;
        astro->sid[0] = ad.sidtime;
        astro->dra[0] = ad.dra;
        astro->rsum[0] = ad.rsum;


    } else if (astro->jd != julianDay)
    {
        
        computeAstroDay(julianDay-1, &ad);
        astro->ra[0] = ad.ra;
        astro->dec[0] = ad.dec;
        astro->sid[0] = ad.sidtime;
        astro->dra[0] = ad.dra;
        astro->rsum[0] = ad.rsum;
        computeAstroDay(julianDay, &ad);
        astro->ra[1] = ad.ra;
        astro->dec[1] = ad.dec;
        astro->sid[1] = ad.sidtime;
        astro->dra[1] = ad.dra;
        astro->rsum[1] = ad.rsum;
        computeAstroDay(julianDay+1, &ad);
        astro->ra[2] = ad.ra;
        astro->dec[2] = ad.dec;
        astro->sid[2] = ad.sidtime;
        astro->dra[2] = ad.dra;
        astro->rsum[2] = ad.rsum;
    
    }

    astro->jd = julianDay;
    computeTopAstro(loc, astro, topAstro);

}


void computeAstroDay(double JD, AstroDay* astroday)
{

    int i =0;
    double R, Gg, G;

    double tL, L;
    double tB, B;

    double X0, X1, X2, X3, X4;

    double U, E0, E, lamda, V0, V;

    double RAn, RAd, RA, DEC;

    double B0sum=0, B1sum=0;
    double R0sum=0, R1sum=0, R2sum=0, R3sum=0, R4sum=0;
    double L0sum=0, L1sum=0, L2sum=0, L3sum=0, L4sum=0, L5sum=0;

    double xsum=0, psi=0, epsilon=0;
    double deltaPsi, deltaEps;

    double JC = (JD - 2451545)/36525.0;                                                 
    double JM = JC/10.0 ; 
    double JM2 = pow (JM, 2);
    double JM3 = pow (JM, 3);
    double JM4 = pow (JM, 4);
    double JM5 = pow (JM, 5);



    for(i=0; i < 64; i++)
        L0sum += L0[i][0] * cos(L0[i][1] + L0[i][2] * JM);
    for(i=0; i < 34; i++)
        L1sum += L1[i][0] * cos(L1[i][1] + L1[i][2] * JM);
    for(i=0; i < 20; i++)
        L2sum += L2[i][0] * cos(L2[i][1] + L2[i][2] * JM);
    for(i=0; i < 7; i++)
        L3sum += L3[i][0] * cos(L3[i][1] + L3[i][2] * JM);
    for(i=0; i < 3; i++)
        L4sum += L4[i][0] * cos(L4[i][1] + L4[i][2] * JM);
    L5sum = L5[0][0] * cos(L5[0][1] + L5[0][2] * JM);

    
    tL = (L0sum + (L1sum * JM) + (L2sum * JM2) 
          + (L3sum * JM3) + (L4sum * JM4) 
          + (L5sum * JM5)) / pow (10, 8);
 
    L = limitAngle(RAD_TO_DEG(tL));

    for(i=0; i<5; i++)
        B0sum += B0[i][0] * cos(B0[i][1] + B0[i][2] * JM);
    for(i=0; i<2; i++)
        B1sum += B1[i][0] * cos(B1[i][1] + B1[i][2] * JM);


    tB= (B0sum + (B1sum * JM)) / pow (10, 8);
    B = RAD_TO_DEG(tB);

    
    for(i=0; i < 40; i++)
        R0sum += R0[i][0] * cos(R0[i][1] + R0[i][2] * JM);
    for(i=0; i < 10; i++)
        R1sum += R1[i][0] * cos(R1[i][1] + R1[i][2] * JM);
    for(i=0; i < 6; i++)
        R2sum += R2[i][0] * cos(R2[i][1] + R2[i][2] * JM);
    for(i=0; i < 2; i++)
        R3sum += R3[i][0] * cos(R3[i][1] + R3[i][2] * JM);
    R4sum = R4[i][0] * cos(R4[i][1] + R4[i][2] * JM);

    R = (R0sum + (R1sum * JM) + (R2sum * JM2)
         + (R3sum * JM3) + (R4sum * JM4)) / pow (10, 8);

    G = limitAngle((L + 180));
    Gg = -B;

    X0 = 297.85036 + (445267.111480 * JC) -  (0.0019142 * pow (JC, 2)) +
        pow (JC, 3)/189474.0;
    X1 = 357.52772 + (35999.050340 * JC) -  (0.0001603 * pow (JC, 2)) -  
        pow (JC, 3)/300000.0;
    X2 = 134.96298 + (477198.867398 * JC) +  (0.0086972 * pow (JC, 2)) +  
        pow (JC, 3)/56250.0;
    X3 = 93.27191 + (483202.017538 * JC) -  ( 0.0036825 * pow (JC, 2)) +  
        pow (JC, 3)/327270.0;
    X4 = 125.04452 - (1934.136261 * JC) + (0.0020708 * pow (JC, 2)) +  
        pow (JC, 3)/450000.0;


    for (i=0; i<63; i++) {
        xsum += X0*SINCOEFF[i][0];
        xsum += X1*SINCOEFF[i][1];
        xsum += X2*SINCOEFF[i][2];
        xsum += X3*SINCOEFF[i][3];
        xsum += X4*SINCOEFF[i][4];
        psi     += (PE[i][0] + JC*PE[i][1])*sin(DEG_TO_RAD(xsum));
        epsilon += (PE[i][2] + JC*PE[i][3])*cos(DEG_TO_RAD(xsum));
        xsum=0;
    }

    deltaPsi = psi/36000000.0;
    deltaEps = epsilon/36000000.0;

    U = JM/10.0;
    E0 = 84381.448 - 4680.93 * U - 1.55 * pow(U,2) + 1999.25 * pow(U,3) 
        - 51.38 * pow(U,4)  - 249.67 * pow(U,5) - 39.05 * pow(U,6) + 7.12 
        * pow(U,7) + 27.87 * pow(U,8) + 5.79 * pow(U,9) + 2.45 * pow(U,10);
    E = E0/3600.0 + deltaEps;
    lamda = G + deltaPsi + (-20.4898/(3600.0 * R));

    V0 = 280.46061837 + 360.98564736629 * ( JD - 2451545) +  
        0.000387933 * pow(JC,2) - pow(JC,3)/ 38710000.0;
    V = limitAngle(V0) + deltaPsi * cos(DEG_TO_RAD(E));

    RAn = sin(DEG_TO_RAD(lamda)) * cos(DEG_TO_RAD(E)) - 
        tan(DEG_TO_RAD(Gg)) * sin(DEG_TO_RAD(E));
    RAd = cos(DEG_TO_RAD(lamda));
    RA = limitAngle(RAD_TO_DEG(atan2(RAn,RAd)));

    DEC = asin( sin(DEG_TO_RAD(Gg)) * cos(DEG_TO_RAD(E)) + 
                            cos(DEG_TO_RAD(Gg)) * sin(DEG_TO_RAD(E)) * 
                            sin(DEG_TO_RAD(lamda)));

    astroday->ra = RA;
    astroday->dec = DEC;
    astroday->sidtime = V;
    astroday->dra = 0;
    astroday->rsum = R;

}

void computeTopAstro(const Location* loc, const Astro* astro, Astro* topAstro)
{
    int i;
    double lHour, SP;
    double tU, tCos, tSin, tRA0 ,tRA ,tDEC;

    for (i=0; i<3; i++)
    {
        lHour = limitAngle(astro->sid[i] + loc->degreeLong - astro->ra[i]);

        SP = 8.794/(3600 * astro->rsum[i]);

        tU = atan (0.99664719 * tan(DEG_TO_RAD(loc->degreeLat)));

        tCos = cos(tU) + ( (loc->seaLevel)/6378140.0) * 
            cos(DEG_TO_RAD(loc->degreeLat));

        tSin = 0.99664719 * sin(tU) + ( loc->seaLevel/6378140.0) * 
            sin(DEG_TO_RAD(loc->degreeLat));

        tRA0 = (((-tCos) * sin(DEG_TO_RAD(SP)) * sin(DEG_TO_RAD(lHour))) 
                / (cos(astro->dec[i]) - tCos * sin(DEG_TO_RAD(SP)) * cos(DEG_TO_RAD(lHour))));


        tRA = astro->ra[i] +  RAD_TO_DEG(tRA0);


        tDEC = RAD_TO_DEG(atan2((sin(astro->dec[i]) - tSin * sin(DEG_TO_RAD(SP))) * cos(tRA0), 
                                cos(astro->dec[i]) - tCos * sin(DEG_TO_RAD(SP)) * 
                                cos(DEG_TO_RAD(lHour))));

        topAstro->ra[i] = tRA;
        topAstro->dec[i] = tDEC;
        topAstro->sid[i] = astro->sid[i];
        topAstro->dra[i] = tRA0;
        topAstro->rsum[i] = astro->rsum[i];

    }

}

double limitAngle(double L)
{
    double F;
    L /= 360.0;
    F = L - floor(L);
    if (F > 0)
        return 360 * F;
    else if (F < 0)
        return 360 - 360 * F;
    else return L; 
}


double limitAngle180(double L)
{
    double F;
    L /= 180.0;
    F = L - floor(L);
    if (F > 0)
        return 180 * F;
    else if (F < 0)
        return 180 - 180 * F;
    else return L;
}

double limitAngle111(double L)
{
    double F;
    F = L - floor(L);
    if (F < 0)
        return F += 1;
    return F;
}

double limitAngle180between(double L)
{
    double F;
    L /= 360.0;
    F = (L - floor(L)) * 360.0;
    if  (F < -180)
        F += 360;
    else if  (F > 180)
        F -= 360;
    return F;
}
    
